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ABSTRACT 

Context. As the coronal magnetic field can usually not be measured directly, it has to be extrapolated from photospheric measurements into 
the corona. 

Aims. We test the quality of a non-linear force-free coronal magnetic field extrapolation code with the help of a known analytical solution. 
Methods. The non-linear force-free equations are numerically solved with the help of an optimization principle. The method minimizes an 
integral over the force-free and solenoidal condition. As boundary condition we use either the magnetic field components on all six sides of 
the computational box in Case I or only on the bottom boundary in Case II. We check the quality of the reconstruction by computing how well 
force-freeness and divergence-freeness are fulfilled and by comparing the numerical solution with the analytical solution. The comparison is 
done with magnetic field line plots and several quantitative measures, like the vector correlation, Cauchy Schwarz, normalized vector error, 
mean vector error and magnetic energy. 

Results. For Case I the reconstructed magnetic field shows good agreement with the original magnetic field topology, whereas in Case II there 
are considerable deviations from the exact solution. This is corroborated by the quantitative measures, which are significantly better for Case I. 
£j Conclusions. Despite the strong nonlinearity of the considered force-free equilibrium, the optimization method of extrapolation is able to 
. . i reconstruct it; however, the quality of reconstruction depends significantly on the consistency of the input data, which is given only if the 
known solution is provided also at the lateral and top boundaries, and on the presence or absence of flux concentrations near the boundaries of 
the magnetogram. 
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1. Introduction 

Several methods have been proposed to compute the non- 
linear force-free coronal magnetic field in active regions from 
measurements of the photospheric mag netic field vector , e.g. 
the direct upward integra tion method dWu et all Il990h. the 



Grad-Rubin method (e.g. ISa kurai 
Reenier et all l2002t iBlevbel et al 



£3 



2002; 



lax method ( Roumeliotisl 199rj; 



Valori et al 



Wiegelmann et al 1 l2006h 



Amari et al. 
IWheatlandl 



1999; 
2004; 



Amar i et al Il2006t llnhester & Wiegelmann|[2"o 06). the Green's 
function like methods jYan & Sakurai 2000), the stress and re 



20051) and the op- 



timization method (IWheatland et al ■I l2000i: IWiegelmannll2004t 
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A standard test for force-free extrapolation methods is the 
application of the codes to known analytical or numerical non- 
linear force-free equilibria. Due to the general mathematical 
difficulty of the problem, only very few such solutions are 
known. The comparison of numerical nonlinear force-free ex- 
trapolation cod es with the c l ass of nonlinear force-free equi- 
libria found by iLow & Loul (119901) (LL from here on) has in 
the past ten years emerged as a certain standard test. The LL 
equilibria are a class of axisymmetric equilibria which are sep- 
arable in spherical coordinates. They are self-similar in the ra- 
dial coordinate, and the polar angle dependence is determined 
from a non-linear eigenvalue equation. The symmetry is bro- 
ken by cutting out a rectangular chunk of the solution by using 
a Cartesian coordinate system which is shifted and rotated with 
respect to the original coordinate system in which the LL equi- 
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libria are calculated. The parameters of the LL solutions and 
the parameters of the new Cartesian coordinate system allow 
for a large number of different situations which can be used for 
tests. 

Nevertheless, it is highly desirable to also use tests differ - 
ent from LL. A first step has been made by lValori et al.l (120051) . 
who used an equilibrium taken from a numerical investigation 
of a twisted flux rope. Valori and Kliem (2005, private com- 
munication) have further tested their nonlinear force-free ex- 
trapolation code with configurations containing a flux rope by 
considering the approximate analytical fo rce-free equilibrium 
established by (iTitov & Demoulinl 1 19991 henceforth TD). In 
a recent review talk on nonlinear force-free extrapolation, T. 
Neukirch also strongly recommended the use of other equilib- 
ria, in particular the TD-equilibrium, to test the quality of non- 
linear force-free extrapolation codes. The TD equilibrium has a 
much more concentrated current density than LL, and the mag- 
netogram is more structured. Also, the final equilibrium has a 
different topology than the corresponding potential field, which 
is not true for LL. 

The aim of this research note is to study the perfor- 
mance the optimization method in reconstructing the TD equi- 
librium. The optimization code has so far only been tested 
with th e LL solution wi t h rathe r small resolution (40 x 40 x 
20) in I Wheatland et al.l d2000h on 80 x 80 x 40 grids in 
Wiegelmann & Neukirch! ( 12003b and on a 64 3 grid as part 
of a comparison of six different extrapolation methods in 
Schrijver et al.l ( 2006b. The op timization code in the implemen- 
tation of IWiegelmannl d2004b was the fastest-converging and 
best-performing model of the six tested extrapolation codes. 
Here we test this optimization code with the TD equilibrium 
using a high numerical resolution (150 x 250 x 100). 

2. The TD equilibrium 

Like the LL solutions, the TD equilibrium is an axisymmet- 
ric equilibrium and the symmetry is broken by choosing the 
boundaries appropriately. The TD equilibrium has originally 
been constructed as a model of an active region containing a 
current- carrying flux tube. The model has later been used in 
a number of investigations of the initiation of CME eruptions 



Table 1. Parameter set of the TD equilibrium 



a nu mber ot inve stigations ot the initiation ot livib eru ptions 
(e.g.jRoussev et alj2003l:lTorok et allfooH lKliem et all2004 



Torok & Kliemll2005HWilliams et al.ll2005l) 

The axis of symmetry of the TD equilibrium is placed at 
a distance d beneath the lower boundary (the photosphere). 
A line current of strength Iq runs along the line of symme- 
try, which creates a potential m agnetic field with circular field 
lines. Titov & Demoulin ( 1999b then added a toroidal nonlinear 
force-free current around the axis of symmetry with minor ra- 
dius a and major radius R and total current /, assuming a «: R 
(see also Fig. 2 in ITitov & Demoulinldl999h ). Similar to a toko- 
mak, the ring current would not be in equilibrium with the po- 
tential field created by the line current. Therefore, two magnetic 
monopoles (strength q) of opposite polarity are placed on the 
symmetry axis at a distance L from each other, with the force- 
free ring current half-way between them. The monopoles create 
a potential poloidal magnetic field which has field lines lying 
above the force-free ring current and thus holding it down. In 
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this way a stable equilibrium situation can be achieved for cer- 
tain parameter combinations (unstable cases have been investi- 
gated in the studies of CME initiation mentioned above). 

For the present investigation we use the TD equilibrium 
(given by their equations 16-22 and 31) with the set of param- 
eters given in Table Q] This yields a stable equilibrium because 
the twist of the flux rope is clearly sub -critical with respect to 
the kink instability dTorok et alj|2004b and the mo nopoles are 
placed sufficiently distant dTitov & Demoulinl 1999b . In TableQ] 
we have also included the derived parameters of the apex height 
of the ring current above the photosphere {- R - d), the foot- 
point position of the ring current (= y/R 2 - d 2 ), the average 
field line twist in the nonlinear force-free ring current (aver- 
aged over the toroidal cross section), and the normalization of 
the field and of the various length parameters. 

The TD equilibrium has a number of properties which we 
mention here explicitly because they are relevant for magnetic 
field extrapolation methods. The artificial magnetogram of the 
TD equilibrium basically has a large-scale bipolar structure 
(due to the line current), but this structure is modified by the 
monopoles and the ring current to an almost quadrupolar struc- 
ture. Because the line current extends from — oo to oo, the mag- 
netic field on the photosphere does not drop to zero with larger 
distance from the nonlinear force-free region in the direction 
of the line current. This is different from the LL solutions, for 
which the photospheric magnetic field can be made very small 
at the edges of the magnetogram if the computational box is 
chosen large enough. On the other hand, the coronal current 
density (and therefore the region with non- vanishing a on the 
photospheric boundary) is strongly concentrated for the TD 
equilibrium (to the region of the ring current). In the LL case 
a is distributed over extended regions of the lower boundary. 
Finally, as already mentioned above, the TD equilibrium mag- 
netic field has a different topology than the corresponding po- 
tential magnetic field calculated from the B : component of the 
magnetic field on the lower boundary. This is not the case for 
the standard LL cases usually used for testing nonlinear force- 
free extrapolation codes. 

3. Method 

Force-free coronal magnetic fields have to obey the equations 

(VxB)xB = 0, (1) 
V B = 0. (2) 
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We solve these equati ons with the help of an o ptimization prin- 
ci ple, as proposed by [W heatland et al. I (l2000l) and generalized 
by Wiegelmann & rnhesterld2003l) . We define the functional 



Demoulin equilibrium) and b (reconstructed fiel ds). We use ex 



L= j w(x,y,z) [fi~ 2 |(V x B) x B| 2 + |V ■ B| 2 ] d\, 



(3) 



where w(x, y, z) is a weighting function. It is obvious that (for 
w > 0) the force-free Eqs. ([T]|2]i are fulfilled when L equals 
zero. We compute the magnetic field in a box with nx = 150, 
ny = 250 and nz = 100 points. The numerical method works 
as follows. As an initial configuration we compute a potential 
magnetic field in the whole bo x with the help of a Green's 
function method as described in lAlvl (11989b . The next step is 
slightly different for the two cases. 

- Case I: We impose the exact TD magnetic field vector on 
all six boundaries of the computational box. The weighting 
function is w — 1 in the entire box. 

- Case II: We impose the exact TD magnetic field vector 
only on the bottom boundary (photosphere). The boundary 
conditions on the lateral and top boundaries of the com- 
putational box are given by the initial potential field. The 
weighting function is w = 1 in the center 118x218x84 
region and drops to with a cosine profile in a 16-pixel 
boundary layer towards the lateral and top boundaries of 
the computational box. 

In both cases we iterate for the magnetic field inside the compu- 
tational box by minimizing Eq. (0. The program is written in 
C and has been parallelized with OpenMP. The computations 
have been done on 8 Procs. The detai ls of the current im ple- 
mentation of our code are described in lWiegelmannl (I2004I) . 



4. Figures of merit 

To assess the quality of our reconstruction we introduce sev- 
eral figures of merit in Table [2] The first column of the ta- 
ble gives the name of model, column two contains the inte- 
gral L — L\ + L2 defined in Eq. [3] the third column L\ shows 
how good the force-free condition (first term of the integral in 
Eq. [3} is fulfilled and the forth column Ln corresponds to the 
solenoidal condition (second term of the integral in Eq.[3). L, 
L\ and L2 are calculated with the magnetic field normalized 
by the average magnetic field strength on the photosphere and 
the length scale normalized by the average box length. We em- 
phasize, however, that we restrict the integration domain to the 
inner part of the computational box, excluding the boundary 
layer in which w(x, y,z) + 1 . Thus strictly speaking the diag- 
nostic quantities differ from the integrals given in Eq. [3] We 
therefore call them Li nner , L\ j nner and dinner- Columns five to 
eleven contain quantitative measures on how good the recon- 
structed solution agrees with the original Titov & Demoulin 
equilibrium. These fi gures of merit have been introduced by 
Schrijver et al. (2006) to compare the results of six different ex- 
trapolation codes (including the code used here) with the Low 
& Lou solution. In a coordinated studv lAmari et al. (2006) used 
the same figures of merit to compare the reconstruction of two 
Grad-Rubin codes also with the Low & Lou solution. These fig- 
ures quantify the agreement between vector fields B (Titov & 



actly t he same definitions as given in Sect. 4 of ISchriiver et al 



(2006), in order to allow also an comparison of how well our 
co de reconstructs the Lo w & Lou solution (see Tables I and II 
in ISchriiver et al. and the Titov & Demoulin solution 

(Table|2]here). 

Column five contains the L M norm of the divergence of the 
magnetic field 



V B |L=sup|V-B| 

xeV 



(4) 



Column six contains the L m norm of the Lorentz force of the 
magnetic field 



I i xB II 

00 — 

sup |j x B| 

xeV 

Column seven contains the vector correlation 

\l/2 



C V ec =2 Brbi/ Xi |Bi|2 Z |bi|2 



V i 



column eight the Cauchy-Schwarz inequality 
Bi bi 



Ccs 



1 y Bi bi 

N 4-* |Bi||bi| 



(5) 



(6) 



(7) 



where N is the number of vectors in the field. The normalized 
vector error is defined as 



£ N = Jjbi-Bi|/^|Bi|, 

i i 

and the mean vector error as 



Em - 



N E/IBil 



(8) 



(9) 



Columns nine and ten contain E N = 1 - En and E M = 1 — 
Em, respectively. In column eleven we show the total magnetic 
energy of the reconstructed field normalized with the energy of 
the input field 



Z,-|bif 

2/|Bi| : 



(10) 



The two vector fields agree perfectly if the figures of merit 
(C vec , Ccs, E N , E M , e) are unity. In the twelves column we list 
ep which is the magnetic energy normalized by the energy 
of the corresponding potential field. Column thirteen contains 
the number of iteration steps until convergence and column 
fourtenn the computing time on 8 processors. 

5. Results 

In Fig. [T] we compare the original TD equilibrium with a po- 
tential field extrapolation (Potential), a non-linear force-free 
extrapolation with all six boundaries prescribed (Case I) and 
a non-linear force free extrapolation where only the bottom 
boundary has been prescribed (Case II). The potential field 
does obviously not agree with the original. Case I shows a very 
good agreement with the original. Case II reconstructs the mag- 
netic field approximately. Deviations of the magnetic field lines 
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Table 2. Figures of merit for the comparison of our reconstruction with the original TD equilibrium. See Sec.|4]for details. The 
diagnostic was done on the center 118x218x84 grid and a boundary layer of 16 points towards the lateral and top boundaries 
of the computational box has been omitted. 
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Fig. 1. We compare the magnetic field structure of the original Titov Demoulin equilibrium (Titov Demoulin) with a potential field 
extrapolation (Potential), a non-linear force-free extrapolation with all six boundary surfaces prescribed (Case I) and a non-linear 
force-free extrapolation which uses only the photospheric magnetic field data (Case II). The panels show the center 118x218x84 
region of the original 150 x 250 x 100. The colour coding shows the magnetic field strength on the photosphere. We used the 
same start points for the magnetic field line computation in the positive (yellow) magnetic field region on the photosphere in all 
panels. 



between Case II and the original are visible, in particular for 
field lines emerging in the nonlinear force-free region of the 
bottom boundary. 

The visible inspection of Fig. [TJ is supported by the fig- 
ures of merit in Tab. [2] The figures of merit (columns 7-11) 
are for case I better than 1 % for all figures, mostly even bet- 
ter than 0.1%. The integral force and divergence free condi- 



tions (columns 1-3) are even one order of magnitude better than 
the discretization error of the exact solution. The L x norms 
(columns 4 and 5) are of the same order of magnitude as the 
the discretization error. 

The reason is that the optimization code minimizes L with 
respect to the numerical grid. For case II, on the contrary, L 
is one order of magnitude higher than the discretization error 



Wiegelmann et al.: Testing magnetic field extrapolations (RN) 



5 



of the exact solution. The reason is that for case II an inconsis- 
tency exists between the bottom boundary and the other bound- 
aries, where the initial potential field is kept. Correspondingly, 
the figures of merit are less good than in case I. 

The lateral and upper boundaries try to feed information 
about the original potential field into the solution, and because 
the potential field has a different topology, the flux rope of the 
TD equilibrium is not correctly reproduced in Case II. Fig. 2 
(online material) shows that the field lines calculated from start 
points lying on symmetric circles at the footpoints of the flux 
rope do outline a single flux rope only in Case I, but outline 
two flux ropes that separate with height in Case II. 

6. Conclusions 

We used the TD equilibrium to test a non-linear force-free op- 
timization code. In case I we prescribed all six boundaries of 
the computational box and in case II only the bottom boundary. 
For case I we get an almost perfect agreement with the origi- 
nal, and for case II also a reasonable agreement. Case I is the 
real performance test of the code, because only if the correct 
boundary conditions on all boundaries are prescribed we have 
exactly the same physical problem and can expect to find the 
exact solution. For the reconstruction of magnetic fields from 
observed vector magnetograms the results of Case II are more 
important, however, because for real active regions the lateral 
boundaries are unknown. Let us also remark that the figures 
of merit here are comparab le with the figu r es of merit in the 
LL case (Table I rows b) in ISchriiver et al. For cases 



II the results for TD are slightly worse for L, L\ , C vec , E N and 
somewhat better for L2, Ccs, E M compared with the LL case. 
The reconstructed magnetic energy in the TD case discussed 
here is not as accurate (only for case II) as for the LL case. One 
has to consider, however, tha t the magnetic energy of the corre- 
sponding LL field (case II in Schriiver et alj|2006l) is only 1.10 
times the potential field, where the TD equilibrium has a mag- 
netic energy which is 2.36 higher than the potential field. The 
magnetic field strength on the photosphere for the LL field be- 
comes very low close to the lateral boundaries, which is not the 
case for TD. Despite these difficulties, our code reconstructed 
the magnetic field topology approximately correct, even if the 
lateral boundaries are unknown. Current vectormagnetograms 
often have a limited field of view and consequently also signif- 
icant magnetic flux close to the lateral boundaries. A compar- 
iso n of direct measurem ents of chromospheric magnetic loops 
by ISolanki et al.1 ( 120031) with a non-linear force-free extrapo- 
lation from a p hotospheric vectormagne togram with a limited 
field of view in lWiegelmann et alJ (120051) showed a reasonable 
agreement, except for loops close to boundaries of the magne- 
togram. 
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Online material 
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Fig. 2. The TD equilibrium has a 
current ring surrounded by potential 
field. Being force-free this means 
that there is a semi-circular flux 
rope in correspondence of the cur- 
rent ring. Here we check how good 
the extrapolated final topology re- 
produces that. We plot two sets of 6 
field lines, each set defining a flux 
tube. The field lines start from 6 
equally spaced points on two cir- 
cles of radius 0.1, each centered 
at one of the intersections of the 
current ring axis with the magne- 
togram (x = z = 0,y = ±y foot ). 
The aim of such a plot is to visu- 
alize the flux rope that sustains the 
current ring in TD. In the original 
TD field the two flux tubes over- 
lap completely: one ends where the 
other starts. Case I confirms that the 
extrapolation using the information 
on all six boundaries recovers the 
flux rope to an excellent degree of 
accuracy. However, when only one 
boundary is used for the extrapola- 
tion the two flux tubes do not coin- 
cide (Case II). 



